% Simulation
clc; clear all; close all;
location = strcat(pwd,'\');
data = strcat(location,'data\');
figout = strcat(location,'figures\');
path(path,strcat(location,'function_calls'))

% Main Calibration
load(strcat(data,'vf_60_150.mat'));

% Parameters
A = 22; % leverage
K = 5/100; % fixed cost
gamma = 12; % cost convexity parameter (note, as inten<1, higher gamma actually means cheaper)
dt = 1/12;
rho = (5/100)*dt; % intertemporal discount rate

alpha = alpha_60_150; % adjusted to monthly (freq of simulation)
sigma = var_60_150;

n_abs = inten_60_150; % define absolute intensity
ind = (p>0.98);
n_abs(ind) = 0;

% for below analysis smooth intensity and convert mesh- to ngrid
q_ngrid=q'; p_ngrid=p'; inten_ngrid=n_abs'; v_ngrid=v_60_150'; pi_ngrid=pi_60_150';

%% Monte Carlo Simulations
T = 20;
N = 1000;

rng('default')
rng(1) % fix seed

% unconditional prob of good, generate alphas of fund
prob_good = 0.175; 
fund_alpha = binornd(ones(N,1),prob_good);
fund_alpha(fund_alpha==0) = -alpha;
fund_alpha(fund_alpha==1) = +alpha;

% variables of interest (randomize priors in loop for each fund)
phat = zeros(T/dt+1,N); % allocator probability
qhat = zeros(T/dt+1,N); % market probability
real_ret = zeros(T/dt,N); % realized returns
intensity = zeros(T/dt,N); % intensity decision
inv = zeros(T/dt+1,N); % investment indicator

% interpolate intensity
inten_func = griddedInterpolant(q_ngrid,p_ngrid,inten_ngrid);

% interpolate investment decision 
ind = double(v_ngrid==pi_ngrid & p_ngrid>0.1 & q_ngrid<0.9 & p_ngrid-q_ngrid>0.25);
inv_func = griddedInterpolant(q_ngrid,p_ngrid,ind);

% simulation
for n=1:1:N
    phat(1,n) = 0.30+0.15*rand; % time zero prior of allocator
    qhat(1,n) = max(0.05+0.15*rand,0); % time zero prior of market
    
    for t=1:1:T/dt
    shock = randn; % same signal
    
    % market (objective q)
    dxhat = fund_alpha(n) + sigma*shock;
    real_ret(t,n) = dxhat - alpha*(2*qhat(t,n)-1);
        
    dWhat = (1/sigma)*(dxhat-(2*alpha*qhat(t,n)-alpha)); % filtered shock
    qhat(t+1,n) = real(qhat(t,n) + qhat(t,n)*(1-qhat(t,n))*(2*alpha/sigma)*dWhat); % probability of high type
    
        % allocator (subjective p)
        if inv(t,n) == 1
            inv(t+1,n) = 1;

            dxbar = fund_alpha(n) + sigma*shock;
            dWbar = (1/sigma)*(dxbar-(2*alpha*phat(t,n)-alpha)); % no intensity now (already invested in)
            phat(t+1,n) = real(phat(t,n) + phat(t,n)*(1-phat(t,n))*(2*alpha/sigma)*dWbar);
        else
            intensity(t,n) = inten_func(qhat(t,n),phat(t,n)); % equilibrium intensity

            dxbar = fund_alpha(n) + (sigma/sqrt(1+intensity(t,n)))*shock;
            dWbar = (sqrt(1+intensity(t,n))/sigma)*(dxbar-(2*alpha*phat(t,n)-alpha));
            phat(t+1,n) = real(phat(t,n) + phat(t,n)*(1-phat(t,n))*(2*alpha/sigma)*sqrt(1+intensity(t,n))*dWbar);

            inv(t+1,n) = inv_func(qhat(t+1,n),phat(t+1,n));
            if inv(t+1,n)>0.95
                inv(t+1,n)=1;
            else
                inv(t+1,n)=0;
            end
        end    
    end
    
    string = "Fund number" + " " + num2str(n) + " " + "complete!";
    disp(string)
end

% Ever invest indicator
einv = repmat(max(inv),size(inv,1),1);

% remove first few periods (otherwise it dominates in heatmap as concentrated around prior)
cutoff = 2;
phat = phat(cutoff:end,:);
qhat = qhat(cutoff:end,:);
einv = einv(cutoff:end,:);
inv = inv(cutoff:end,:);

% Save simulation output
S.('einvhat') = einv;
S.('phat') = phat;
S.('qhat') = qhat;
S.('inv') = inv;
S.('intenhat') = intensity;

v = v_60_150;
pi = pi_60_150;
inten = inten_60_150;
S.('v') = v; % keep all variables below in mesh grid format
S.('pi') = pi; 
S.('inten') = inten;
S.('p') = p;
S.('q') = q;
S.('fund_alpha') = fund_alpha;
S.('real_ret') = real_ret;

save(strcat(data,'heatmap_data.mat'),'-struct','S');

%% Calibration Table A1 & other simulation summary stats
load(strcat(data,'heatmap_data.mat')); % simulation data

ind = (fund_alpha>0); % good fund
frac_inv = max(inv);

frac_good = sum(frac_inv(ind))/sum(ind); % fraction of bad funds invested in
frac_bad = sum(frac_inv(~ind))/sum(~ind); % fraction of good funds invested in
frac_inv = sum(frac_inv)/N; % fraction of all funds invested in

inv_stats_all = inv(:,(max(einv)==1)); % pull columns of invested funds
inv_stats_all = inv_stats_all(2:end,:)-inv_stats_all(1:end-1,:); % date of investment
inv_stats_all = [zeros(1,size(inv_stats_all,2)); inv_stats_all];

[rows,~,~] = find(inv_stats_all==1);
avg_due_diligence = mean(rows); % average month of investment

inv_stats_alpha = fund_alpha((max(einv)==1));
ind = (rows>= 22); % indicator of established managers

fund_alpha_old = inv_stats_alpha(ind); % fraction of old funds that are good
good_old_frac = (fund_alpha_old>0);
good_old_frac = sum(good_old_frac)/sum(ind); 

fund_alpha_young = inv_stats_alpha(~ind); % fraction of young funds that are good
good_young_frac = (fund_alpha_young>0);
good_young_frac = sum(good_young_frac)/sum(~ind);

real_ret_inv = real_ret.*inv;
real_ret_inv = sum(real_ret_inv);

einv_cross = max(einvhat);
ind = (einv_cross==1);
real_ret_inv = mean(real_ret_inv(ind)-K);
twntym_real_ret = A*(real_ret_inv/(240-avg_due_diligence));

%% Time Series Example (one young & one established)
load(strcat(data,'heatmap_data.mat')); % simulation data

einv_cross = max(einvhat);

ind = (einv_cross==1);
phati = phat(:,ind);
qhati = qhat(:,ind);
inv_dat = 2*inv(:,ind)-0.001;

Er = (phati-qhati)*alpha;
chosen_examples = [13 15]; % established / young example
figure_name = {'A1B','A1A'};

for j=1:size(chosen_examples,2)
    i = chosen_examples(j);

    fig = figure('visible','off');
    pl= plot(1:240,Er(:,i),'k',1:240,inv_dat(:,i),'k--');
    pl(1).LineWidth = 2;
    pl(2).LineWidth = 2;
    ylim([0 0.003]); xlim([1 240]);

    ax = gca;
    ax.FontSize = 14;

    ylabel('Monthly Expected Return')
    xlabel('Month from Due-diligence Start')

    set(fig,'PaperSize',fliplr(get(fig,'PaperSize'))) 
    print(fig,'-fillpage','-dpdf','-r500',strcat(figout,'Figure',figure_name{j}))    
end

%% Heat Map
load(strcat(data,'heatmap_data.mat')); % simulation data

% Condition on invested subsample & due-diligence period
ind = (einvhat==1); % invested
phatp = phat(ind);
qhatp = qhat(ind);
invp = inv(ind);

ind = (invp==0); % due-diligence period
phatp = phatp(ind);
qhatp = qhatp(ind);

% define grid
prange = (-0.005:0.01:1.005)';
qrange = (-0.005:0.01:1.005)';

heatmap = nan(size(prange,1)-1,size(qrange,1)-1);
for pind= 1:size(prange,1)-1
    ind = (phatp>=prange(pind) & phatp<prange(pind+1));
    qtemp = qhatp(ind);
    
    for qind=1:size(qrange,1)-1
        ind = (qtemp>=qrange(qind) & qtemp<qrange(qind+1));
        
        heatmap(pind,qind) = size(qtemp(ind),1); % note meshgrid structure
    end
end

[X,Y] = meshgrid(0:0.01:1, 0:0.01:1);
[X2,Y2] = meshgrid(0:0.0001:1, 0:0.0001:1);

fheatmap = interp2(X, Y, heatmap, X2, Y2, 'linear');

inv_grid = double(v==pi & p>0.1 & q<0.9 & p-q>0.25);
finvment = interp2(X, Y, inv_grid, X2, Y2, 'linear');

M = max(max(fheatmap));
bheatmap = fheatmap/M;

fig = figure;
imagesc(fheatmap)
yw = [0.5000         0         0
    0.5417         0         0
    0.5833         0         0
    0.6250         0         0
    0.6667         0         0
    0.7083         0         0
    0.7500         0         0
    0.7917         0         0
    0.8333         0         0
    0.8750         0         0
    0.9167         0         0
    0.9583         0         0
    1.0000         0         0
    1.0000    0.0417         0
    1.0000    0.0833         0
    1.0000    0.1250         0
    1.0000    0.1667         0
    1.0000    0.2083         0
    1.0000    0.2500         0
    1.0000    0.2917         0
    1.0000    0.3333         0
    1.0000    0.3750         0
    1.0000    0.4167         0
    1.0000    0.4583         0
    1.0000    0.5000         0
    1.0000    0.5417         0
    1.0000    0.5833         0
    1.0000    0.6250         0
    1.0000    0.6667         0
    1.0000    0.7083         0
    1.0000    0.7500         0
    1.0000    0.7917         0
    1.0000    0.8333         0
    1.0000    0.8750         0
    1.0000    0.9167         0
    1.0000    0.9583         0
    1.0000    1.0000         0
    1.0000    1.0000    0.0625
    1.0000    1.0000    0.1250
    1.0000    1.0000    0.1875
    1.0000    1.0000    0.2500
    1.0000    1.0000    0.3125
    1.0000    1.0000    0.3750
    1.0000    1.0000    0.4375
    1.0000    1.0000    0.5000
    1.0000    1.0000    0.5625
    1.0000    1.0000    0.6250
    1.0000    1.0000    0.6875
    1.0000    1.0000    0.7500
    1.0000    1.0000    0.8125
    1.0000    1.0000    0.8750
    1.0000    1.0000    0.9375
    1.0000    1.0000    1.0000];

colormap(flipud(gray)) % use "yw" RGB if you want only yellow-to-white
set(gca,'YDir','normal','XDir','normal')
ylim([0 10001]); xlim([0 10001]);
hold on;

BW = imbinarize(bheatmap); % establish boundaries of active region
B = bwboundaries(BW);
    
inv_boundary = [1,1];
[~,h] = contour(finvment,inv_boundary,'k-.');
h.LineWidth = 4;
set(gca,'visible','on')
xlabel('Market Assess. (P(High | Public))')
ylabel('Allocator''s Assess. (P(High | Public + Private))')

ax = gca;
ax.FontSize = 14;

cm_ticks=0:0.1:1;
px_ticks_x=linspace(0,10001,numel(cm_ticks));
px_ticks_y=linspace(0,10001,numel(cm_ticks));
ticklabels=cellfun(@(v) sprintf('%.1f',v),num2cell(cm_ticks),...
    'UniformOutput',false);
ticklabels{1} = '0.0';

set(gca,'XTick',px_ticks_x)
set(gca,'XTickLabels',ticklabels)
set(gca,'YTick',px_ticks_y)
set(gca,'YTickLabels',ticklabels)

set(fig,'PaperSize',fliplr(get(fig,'PaperSize'))) 
print(fig,'-fillpage','-dpdf','-r500',strcat(figout,'Figure3C'))
hold off;

%% E[r] graph as a function of intensity
load(strcat(data,'heatmap_data.mat')); % simulation data

Er = -K+A*(2*(p-q)*alpha);
y = Er(:);
x = inten(:);

% Condition on simulation output
ind = (inv==0); % due-diligence period
sim_data_all = [round(phat(ind),2) round(qhat(ind),2)];
sim_data = unique(sim_data_all,'rows');

grid = [p(:) q(:)];
ind = ismember(grid,sim_data,'rows'); % state space that is active in simulated data

x_active = x(ind);
y_active = y(ind);
grid_active = grid(ind,:);

count = zeros(size(grid_active,1),1);
for i=1:size(grid_active,1) % 
    count(i) = sum(sim_data_all(:,1)==grid_active(i,1) & sim_data_all(:,2)==grid_active(i,2));
end

cutoff = 0.175; % intensity cutoff at low intensity very non-linear relationship
ind = (x_active>cutoff); % define dd region

x_active = x_active(ind);
y_active = y_active(ind);
count = count(ind);

% cutoff some wierd outlying points (due to numerical issues at p=1 boundary)
m = (-0.0044+0.002)/(0.4838-0.4877); % this is for lower right hand side
b = -0.0044-m*0.4838;
temp = y_active-m*x_active-b; % define plane below which you keep points
ind = (temp>0);
%x_active = x_active(ind);
%y_active = y_active(ind);

% Create bins of equal sizes
temp = sort(x_active); 
edges = temp(1:round(length(temp)/25):end);
edges(end) = temp(end);

avg_ret = zeros(size(edges,2)-1,1); % depository for weighted return & intensity in bin
avg_scale = zeros(size(edges,2)-1,1);

fig = figure;
for i=2:size(edges,1)
    ind = (x_active>edges(i-1) & x_active<=edges(i));
    ret = y_active(ind);
    scale = x_active(ind);
    
    color_spread = count(ind)./max(count(ind));
    %{
    for j=1:size(ret,1)
        s = scatter(scale(j),ret(j),'o','filled');
        s.SizeData = 15;
        s.MarkerEdgeColor = (1-color_spread(j))*[1 1 1];
        s.MarkerFaceColor = (1-color_spread(j))*[1 1 1];
        s.MarkerEdgeAlpha = 0.3;
        s.MarkerFaceAlpha = 0.3;
        hold on;
    end
    %}
    w = color_spread./sum(color_spread);
    avg_ret(i-1) = sum(w.*ret);
    avg_scale(i-1) = sum(w.*scale);
end

P0 = [-0.045 0.04 0.43 0.375];
lb = [-0.060 0.02 0.02 0.2];
ub = [0 0.5 0.8 0.5];

plusfun = @(x) max(x,0);
model = @(P,x) P(1) + P(2)*x + P(3)*plusfun(x-P(4));
Pfit1 = lsqcurvefit(model,P0,x_active,y_active,lb,ub);

scalefit = 0.22:0.001:0.49;
retfit = model(Pfit1,scalefit);

lgray = [0.3810    0.3810    0.3810
    0.3968    0.3968    0.3968
    0.4127    0.4127    0.4127
    0.4286    0.4286    0.4286
    0.4444    0.4444    0.4444
    0.4603    0.4603    0.4603
    0.4762    0.4762    0.4762
    0.4921    0.4921    0.4921
    0.5079    0.5079    0.5079
    0.5238    0.5238    0.5238
    0.5397    0.5397    0.5397
    0.5556    0.5556    0.5556
    0.5714    0.5714    0.5714
    0.5873    0.5873    0.5873
    0.6032    0.6032    0.6032
    0.6190    0.6190    0.6190
    0.6349    0.6349    0.6349
    0.6508    0.6508    0.6508
    0.6667    0.6667    0.6667
    0.6825    0.6825    0.6825
    0.6984    0.6984    0.6984
    0.7143    0.7143    0.7143
    0.7302    0.7302    0.7302
    0.7460    0.7460    0.7460
    0.7619    0.7619    0.7619
    0.7778    0.7778    0.7778
    0.7937    0.7937    0.7937
    0.8095    0.8095    0.8095
    0.8254    0.8254    0.8254
    0.8413    0.8413    0.8413
    0.8571    0.8571    0.8571
    0.8730    0.8730    0.8730
    0.8889    0.8889    0.8889
    0.9048    0.9048    0.9048
    0.9206    0.9206    0.9206
    0.9365    0.9365    0.9365
    0.9524    0.9524    0.9524
    0.9683    0.9683    0.9683
    0.9841    0.9841    0.9841
    1.0000    1.0000    1.0000];

scattercloud(x_active,y_active,25,1,'k',flipud(lgray))
set(gca,'Layer','top') % move image to back, axis to front
x_lb = min(x_active);
x_ub = max(x_active);
y_lb = min(y_active);
y_ub = max(y_active);

hold on;
plot(scalefit,retfit,'k','LineWidth',2)

% phat-qhat, expected returns on investment date
inv_stats_all = inv(2:end,:)-inv(1:end-1,:);
ind = (inv_stats_all==1);

phatt = phat(1:end-1,:);
qhatt = qhat(1:end-1,:);

inv_stats_all = [round(phatt(ind),2) round(qhatt(ind),2)]; % intensity month before investment
inv_stats = unique(inv_stats_all,'rows');

ind = ismember(grid,inv_stats,'rows'); % state space points

x_active = x(ind);
y_active = y(ind);
grid_active = grid(ind,:);

% cutoff some wierd outlying points (due to numerical issues at p=1 boundary)
ind = (x_active>0.1 & x_active<0.42 & y_active>0.01);
x_active = x_active(~ind);
y_active = y_active(~ind);

[count,edges] = histcounts(y_active); % fewer observations so use bins to show which lvl of Er most active

color_spread = count./max(count);
for i=2:size(edges,2)
    ind = (y_active>edges(i-1) & y_active<=edges(i));
    ret = y_active(ind);
    scale = x_active(ind);

    s = scatter(scale,ret,'^','filled');
    s.SizeData = 25;
    s.MarkerEdgeColor = [0 0 0];
    s.MarkerFaceColor = [0 0 0];
    s.MarkerEdgeAlpha = 0.5*color_spread(i-1);
    s.MarkerFaceAlpha = 0.5*color_spread(i-1);
    hold on;
end

ylim([y_lb y_ub]); xlim([x_lb x_ub]);
yticks([-0.07 -0.05 -0.03 -0.01 0.01 0.03 0.05])
xticks([0.2 0.30 0.40 0.50])
xlabel('Intensity')
ylabel('Expected Rate of Return | Investment')
ytickformat('%.2f')
xtickformat('%.2f')
box on;

ax = gca;
ax.XColor = [0 0 0];
ax.YColor = [0 0 0];
ax.FontSize = 16;
ax.LineWidth = 1.5;
set(fig,'PaperSize',fliplr(get(fig,'PaperSize'))) 
print(fig,'-fillpage','-dpdf','-r500',strcat(figout,'Figure2'))
hold off;